#include "BesselFunctions.h"
#include <math.h>
#include "struve0.h"
#include "struve1.h"
#include "kernel/feInterpolation.h"


BesselFunctions::BesselFunctions():  feSingleton<BesselFunctions>()
{
	LastX = -1;
	LastL1 = 0;
	LastL2 = 0;
	Init();
}


//-----------------------------------------------------------------------------------------------------------------------------
void BesselFunctions::Init()
{
	double x = 0;
	double step = J0_MAX/J0_SIZE;
	for(int i=0; i<J0_SIZE; i++)
	{
		J0func[i] = _j0(i*step);
		J1func[i] = _j1(i*step);
		Arg[i] = i*step;
	}
}


//-----------------------------------------------------------------------------------------------------------------------------
double BesselFunctions::J0(double iX)
{
	double step = J0_MAX/J0_SIZE;

	int n = int(iX/step + 0.5);


	//get from table
	if(n > J0_SIZE) return 0;
	if(n == 0) n = 1;
	if(n == J0_SIZE-1) n--;

	double l1 = LastL1;
	double l2 = LastL2;
	if(iX != LastX)
	{
		LastX = iX;
		LastL1 = l1 = feInterpolation::LinearBasis1(Arg[n], Arg[n+1], iX);
		LastL2 = l2 = feInterpolation::LinearBasis2(Arg[n], Arg[n+1], iX);
	}
	return J0func[n]*l1 + J0func[n+1]*l2;
}


//-----------------------------------------------------------------------------------------------------------------------------
double BesselFunctions::J1(double iX)
{
	double step = J0_MAX/J0_SIZE;

	int n = int(iX/step + 0.5);
	
	if(n > J0_SIZE) return 0;

	if(n == 0) n = 1;
	if(n == J0_SIZE-1) n--;

	double l1 = LastL1;
	double l2 = LastL2;
	if(iX != LastX)
	{
		LastX = iX;
		LastL1 = l1 = feInterpolation::LinearBasis1(Arg[n], Arg[n+1], iX);
		LastL2 = l2 = feInterpolation::LinearBasis2(Arg[n], Arg[n+1], iX);
	}

	
	return J1func[n]*l1 + J1func[n+1]*l2;
}


//-----------------------------------------------------------------------------------------------------------------------------
double BesselFunctions::StruveH0(double iX)
{
	double step = J0_MAX/J0_SIZE;

	int n = int(iX/step + 0.5);


	//get from table
	if(n > J0_SIZE) return 0;
	if(n == 0) n = 1;
	if(n == J0_SIZE-1) n--;


	double l1 = LastL1;
	double l2 = LastL2;
	if(iX != LastX)
	{
		LastX = iX;
		LastL1 = l1 = feInterpolation::LinearBasis1(Arg[n], Arg[n+1], iX);
		LastL2 = l2 = feInterpolation::LinearBasis2(Arg[n], Arg[n+1], iX);
	}

	return ::StruveH0[n]*l1 + ::StruveH0[n+1]*l2;
}


//-----------------------------------------------------------------------------------------------------------------------------
double BesselFunctions::StruveH1(double iX)
{
	double step = J0_MAX/J0_SIZE;

	int n = int(iX/step + 0.5);


	//get from table
	if(n > J0_SIZE) return 0;
	if(n == 0) n = 1;
	if(n == J0_SIZE-1) n--;

	double l1 = LastL1;
	double l2 = LastL2;
	if(iX != LastX)
	{
		LastX = iX;
		LastL1 = l1 = feInterpolation::LinearBasis1(Arg[n], Arg[n+1], iX);
		LastL2 = l2 = feInterpolation::LinearBasis2(Arg[n], Arg[n+1], iX);
	}

	return ::StruveH1[n]*l1 + ::StruveH1[n+1]*l2;
}


